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Abstract 

The energy and structure of dilute hard- and soft-sphere Bose gases are systematically studied 
in the framework of several many-body approaches, as the variational correlated theory, the Bo- 
goliubov model and the uniform limit approximation, valid in the weak interaction regime. When 
possible, the results are compared with the exact diffusion Monte Carlo ones. A Jastrow type 
correlation provides a good description of the systems, both hard- and soft-spheres, if the hyper- 
netted chain energy functional is freely minimized and the resulting Euler equation is solved. The 
study of the soft-spheres potentials confirms the appearance of a dependence of the energy on 
the shape of the potential at gas paremeter values of x ~ 0.001. For quantities other than the 
energy, such as the radial distribution functions and the momentum distributions, the dependence 
appears at any value of x. The occurrence of a maximum in the radial distribution function, in the 
momentum distribution and in the excitation spectrum is a natural effect of the correlations when 
x increases. The asymptotic behaviors of the functions characterizing the structure of the systems 
are also investigated. The uniform limit approach results very easy to implement and provides a 
good description of the soft-sphere gas. Its reliability improves when the interaction weakens. 

PACS numbers: 03.75.Hh, 05.30.Jp, 67.40.Db 
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I. INTRODUCTION 



The study of dilute systems has met a surge of renewed interest in the last years, following 
the experimental achievement of Bose-Einstein condensation in low density atomic gases 
confined in harmonic traps. Systems are termed as dilute when the average interparticle 
distance, r av , is much larger than the range of the interaction, tq. The main parameter 
characterizing the interaction in the dilute regime is the s-wave scattering length, a, and 
the diluteness condition can be expressed in terms of the density, p, as x « 1, where 
x = pa 3 is the gas parameter. 

Positive values of the scattering length correspond to an essentially repulsive interaction 
at short distances. In fact, for an infinite repulsive barrier, with no attractive part, a is just 
given by the radius of the barrier itself. The hard sphere (HS) potential is largely adopted to 
study low density (LD) gases with positive scattering lengths because of its formal simplicity. 
However, other potentials, providing the same value of a, are indistinguishable from the HS 
one at very low x (universality property). Additional details of the interaction, giving raise 
to nonuniv ers al effects, become relevant when the density (and the gas parameter) increases, 
allowing for discriminating among different potential shapes. Typical ranges of x where these 
differences begin to show up have been founcU^ to be x > 10~ 3 . 

In the density ranges attained in Bose-Einstein condensation experiments the local gas 
parameter may well exceed this value, exploiting the large variation of the scattering length 
in the vicinity of a Fehsbach resonance^. In order to quantitatively study this regime it is 
compulsory to check the reliability of the theories adopted in the analysis. A first, and much 
needed step in this direction, consists in understanding the properties of the underlying 
homogeneous gas of bosons at T=0 temperature. This is a time honored subject of many- 
body physics, addressed by a variety of methods, such as perturbative expansions, Monte 
Carlo (MC) type samplings, variational theories, and so on. 

Expansions in x may be suited to study dilute systems whose interaction can be safely 
described in terms of the s-wave scattering length. They can be derived in the framework of 
standard perturbation theories, built on the basis constituted by the ground state of a non 
interacting Bose fluid and its excited states, obtained by promoting particles from the zero 
momentum condensed state to the non-zero ones. Infinite sums of ladder diagrams must 
be accomplished for strongly repulsive potentials. This procedure results in the well known 
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Lee and Yang low density expansion for the energy per particle of a homogeneous gas of 
bosons^, 

E ( h 2 \ A r 128 fx] 

iV = (wj 47ra; i 1 + l5-V^- (1) 
Low order perturbative diagrams can be retained in the weakly interacting Bose gas^ char- 
acterized by an interaction having a finite Fourier transform. 

Alternatively, the interaction may be mimicked in the low density region by a delta- 
shaped pseudopotential 6 , related to a by 

V(r) = Anah 2 5{r)/m , (2) 

which reproduces the value of the scattering length in the Born approximation: 

This approach is followed in the Bogoliubov theory 7 , by introducing a model hamiltonian 
containing the pseudopotential. A canonical transformation of the BCS type allows for a 
partial diagonalization of the model hamiltonian, leading to a x-expansion as in ((TJ) for the 

ground-state energy. 

Nonuniversal effects have been studied in effective field theory in Ref.—, and found to 
depend, at the leading orders, on the effective range, r s , and on a three-body contact force 
parameter, c, determined in this reference by comparison with the exact diffusion Monte 
Carlo results of Ref. 1 . 

The bosonic many-body Schroedinger equation may be solved, for any potential, by 
Monte Carlo based methods, as the Green's Function Monte Carlo and the Diffusion Monte 
Carlo (DMC) ones. Both approaches are exact, apart from statistical errors, and provide 
essential benchmarks to test the reliability of other theories, at least in those cases where 
the numerical accuracy allows for an unbiased analysis of the results. 

The variational approach is carried on within the Correlated Basis Functions (CBF) 
perturbation theory, and represents its zeroth order. The correlated ground state wave 
function for N interacting particles is obtained by acting with a many-body correlation 
operator on the non-interacting ground state wave function 9 : 

* (l,2,...,jV) =^(l,2,...,A0$ o (l,2,...,A0. (4) 

The operator 2, . . . , N) is meant to take into account the spatial correlations induced by 
the interaction on the free wave function $o(l> 2, . . . , N) (for homogeneous bosons $o — !)• 
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For instance, in the HS case the correlation operator prevents the distance between any 
pair of particles from being smaller than the core radius, so that the wave function vanishes 
for these configurations. The theory is variational in the sense that T is determined by 
minimizing the ground state energy. A correlated perturbation theory may be constructed 
by applying the correlation operator to the non-interacting excited states. 

The weakly interacting Bose fluid has been recently studied within the variational method 
by means of the Independent Pair Correlations (IPC) approach^. The IPC wave function is 
written as 

* JPC (1, 2, . . . , N) = 1 + X; h(r tJ ) + X - Knj)h(r lm ) + . . . , (5) 

i<3 (i<j)^(l<m) 

where the two-body correlation functions, h(r), always act on non-overlapping, independent 
pairs. The structure of a Bose system described by the IPC wave function has been analyzed 
by means of an expansion in cluster diagrams. Renormalized HyperNetted Chain (HNC) 
integral equations exactly sum all the diagrams and have been used to compute the energy, 
distribution function, momentum distribution and pairing function of the bulk boson system 
of soft spheres. The Bogoliubov transformation itself is known to generate a wave function 
having the IPC structure^, and a careful comparison between the variational and Bogoliubov 
approaches has been performed in Ref.— , as well as with the outcomes of the Diffusion Monte 
Carlo and low density theories, for soft sphere (SS) and gaussian like potentials. 

The IPC wave function fails for a hard-sphere fluid, since all pairs need to be simultane- 
ously correlated. This task may be achieved by a Jastrow correlated wave function^ 

* J (l,2,...,N) = J[[l+h(r ij )], (6) 

i<j 

where the productory runs over all pairs. The Jastrow ansatz is clearly richer than the 
IPC one, however it has the disvantage that the HNC equations can be only approximately 
solved, since a class of diagrams (the elementary diagrams) is not summable in a closed 
wayiS. The approximations in the solution of the Jastrow-HNC equations are not expected 
to be quantitatively relevant in the low density regime, however the computed energy is, in 
principle, no longer a rigorous upperbound to the true one, in contrast with the IPC case. 

In this paper we consider a system of N spinless bosons of mass m in a volume Q, 
described by the hamiltonian 

t2 JV N 

tf = -^Ev? + E^«) (7) 

3 =l l=t<j 
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where V(r) is a two-body, spherically symmetric potential, in the thermodynamic limit ( iV 
and Q — ► oo, keeping the density, p = N/Q, constant). 

Two different representative choices for the potential are studied: the hard-sphere poten- 
tial, 

!oo r < a 
(8) 
r > a , 

where the diameter a of the hard sphere coincides with the s-wave scattering length, and a 
soft sphere potential, 

f > r < R 
V(r) = \ (9) 

[ r > R , 

whose s-wave scattering length is given by 

a = R[l- t&nh(K R)/(K R)} , (10) 

with Kl = V m/h 2 . 

The optimal Jastrow correlated wave function, obtained by the minimization of the en- 
ergy per particle through the solution of the HNC Euler equation^, is used along this work. 
We compute the radial distribution function, g(r), the static structure function, S(k), and 
the momentum distribution, n(k), for HS and SS models with identical scattering lengths, 
to ascertain the dependence on the potential form along the gas parameter. We work in the 
HNC/0 approximation, which amounts to disregard the contribution from the elementary 
diagrams and whose accuracy is tested by comparison with the exact DMC results. More- 
over, the Jastrow theory is compared with the LD expansion and with the Bogoliubov and 
IPC theories for the SS case. In addition, the SS potential is examined in the so called 
uniform limit (UL)-, defined by assuming \g(r) — 1| << 1 for all r. Special emphasis is 
devoted to the analysis of several asymptotic behaviors. This is an interesting issue that 
cannot be fully addressed with DMC methods, since the limited size of the simulation box 
strongly limits the possibility of studying effects related to long range structures. 

The plan of the paper is as follows: the HNC/0 and Euler equations are shortly revisited 
in Section II; the uniform limit for the SS potential and its connection with the Random 
Phase Approximation (RPA) theory are discussed in Section III; Section IV present the 
results for the HS and SS models and summary and conclusions are given in Section V. 
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II. HNC THEORY 



The exact wave function of a homogenous, interacting Bose system can be written as the 
product of up to iV-body correlation factors 9 

N N 

%{l,2,...,N) = l[f 2 {r ij ) J] f3(r k i,r km ,T lm )--- , (11) 

i<j k<l<m 

where = | — Vj | and f n >0 for all n. High density and strongly interacting systems, like 
atomic liquid 4 He, are accurately described by keeping only two- and three-body correlations 
(fn>3 = !)• Moreover, a proper choice of f 2 already contributes by ~ 70% to the 4 He 
energy 10 . So, for weak interactions and/or low density fluids, the simpler Jastrow correlated 
wave function of Eq.© may largely be enough to capture the essential features of the exact 
ground state and to provide a quantitatively correct description. 



A. Radial distribution function and Euler equation 



The optimal Jastrow correlation function is obtained by minimizing, without restrictions, 
the expectation value of the hamiltonian (|7J), giving the ground state energy, 
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This is accomplished by solving the Euler-Lagrange equatio^ 1 



(12) 



0. 



(13) 



The link between the energy and f 2 is provided by the radial distribution function, g(r) 
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as, in fact, 



^ = \p J dv 12 g(r 12 ) 



V(r 12 )-^V 2 ln/ 2 (r 12 ) 
2m 



(14) 



(15) 



The Jackson-Feenberg identity 9 - has been used to derive Eq. ([TK|) . In turn, the radial distri- 
bution function can be computed by solving the HyperNetted Chain equations: 



9(r 12 ) = f 2 {r l2 )e N ^ +E ^ 
N(r 12 ) = p f dr 3 [g(rn) - lMr w ) - 1 - N(r 32 )\ , 



(16) 



where N(r) and E(r) are the sum of the nodal and elementary diagrams, respectively^. 
The function E(r) is an input to the theory, and the solution of the HNC equations depends 
on its choice. In the HNC/0 scheme this function is set to be zero. This seemingly drastic 
approximation is, however, reliable at low densities since the diagrams contributing to E(r), 
being highly interconnected, are relevant mostly in the large density regions. 

By inverting the relations (|16|) it is possible to express the energy as a functional of the 
radial distribution function, and formally rewrite the Euler equation as: 

3?-- 

If this procedure is carried on in r-space, the outcome is an integro-differential equation for 
g(r)2&. Alternatively, the energy can be written in terms of the static structure function, 
S(k), defined as the Fourier transform of the radial distribution function, 

S(k) = 1 + p J dre lk - r {g{r) - 1) . (18) 

Variation of E[S] leads to the equations^: 

S(k) = t{k) , (19) 

\Jt 2 {k) + 2V ph (k)t(k) 

with t(k) = h 2 k 2 /2m and V p h{k) the so called particle-hole interaction. In r-space and in 
HNC/0, 

h 2 i — 2 

V ph (r) = g(r)V(r) + - V yjg(r) + [g(r) - 1] uj(r) , (20) 



where the fc-space induced interaction, Uj(k), is 



urn = — 2 m ^ . (2i) 

Eqs. ()19p . ()20p and ()2ip are a set of nonlinear coupled equations to be solved iteratively. 
Finally, the knowledge of S(k) (or g(r)) allows to find the optimal Jastrow correlation 
function by inversion of the HNC/0 equations. 

In the HS model the Euler formalism simplifies, since V p h{r < a) = —uoj{r) and V p h{r > 

2 



a) = h 2 /m V^JgJr) + [g(r) - l)ui(r) 



B. Momentum distribution 



The one-body density matrix for an homogeneous Bose fluid, 

/ \ / \ A JdT 2 dv 3 ---dr N q (l,2,...,N)q (l',2,...,N) 
Mr l5 = P Mv) = N - f - - - , (22) 
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contains essential information about the depletion of the condensate in interacting systems 
and the consequent finite occupation of single particle states carrying non-zero momentum. 
Its diagonal part coincides with the one-body density, and in an homogeneous system pi(r = 
0) = p. The condensate fraction, n , (i.e. the fractional occupation of the k = momentum 
state) is related to the long range order of the one-body density matrix by: n Q = pi(r — > 
oo)/p. 

The associated momentum distribution, n(k), is obtained through the Fourier transform 
of the one-body density matrix, 

n(k) = (2ir) 3 pn 6(k) + J dr exp (ik • r) [p t (r) - pn ] . (23) 

The momentum distribution is normalized as 

1 = T^vT / dkn(k) , (24) 



(2tt) 3 p 

while the kinetic energy per particle can be obtained by n(k) through 



' ' r dk — n(k) . (25) 



N (2vr) 3 p J 2m 

The HNC/0 scheme has been extended to evaluate the one-body density matrix for a 
Jastrow correlated wave function^ 3 -. As a consequence, one gets 

EM = rioe ^r) (26) 

p 

where the new nodal function, N ww (r), is given by 

N ww (r 12 ) = p J dr 3 [g wd (r 13 ) - l][gdw(r 32 ) - 1 - N dw (r 32 )] . (27) 

In the HNC/0 scheme, the functions g wd and N wd are solutions of the set of coupled equations 

9wd(r) = f(r)e N ^ 
N wd (r l2 ) = pj dr 3 [g wd (r l3 ) - l][g(r 32 ) - 1 - N(r 32 )] , (28) 

with g wd {r n ) = g dw ( r 2i) and N wd (r 12 ) = N dw (r 21 ). 

The condensate fraction, uq, is given in terms of the vertex factors, Rw and R d , as 

n = e 2Rw ~ Rd , (29) 

where 

R w = p J dr [g wd (r) - 1 - N wd (r)} - | J dr [g wd (r) - 1] N wd (r) , (30) 
and R d is obtained by substituting in (j3~Uj) g wd — > g and N wd — > N. 



III. SOFT SPHERES IN THE UNIFORM LIMIT APPROXIMATION 



The formalism presented in the previous section is independent on the shape of the po- 
tential. Hence it is equally well suited to analyze both HS and SS systems. However, 
everywhere bounded potentials, as the SS one, allow for alternative calculations comple- 
menting the HNC results. A simple estimate of the energy, in this case, is provided by first 
order perturbation theory, which yields an upper-bound to the exact ground state energy: 

N N 2 y 3 2 v ' ' y ' 

where V is the Fourier transform of the potential and <pQ = 1/VL N / 2 is the wave function of 
the corresponding free system, with all particles occupying the zero momentum state. The 
second order perturbative correction to the energy is also easily obtained as 

E 2 (p) 1 f dq | V(q) | 2 



N 2pJ (2tt) 3 h 2 q 2 /m ' (32) 

however the resulting total energy E(p) = Ei(p) + E 2 (p) is no longer a bound to the exact 
one. 

If the interaction is finite the uniform limit approximation^ may also give an accurate 
description of a dilute system. In this regime, correlations are assumed to be weak and the 
ground state is only slightly affected by them. This condition reads, in terms of the radial 
distribution function, as \g{r) — 1| << 1, as already mentioned in the Introduction. In the 
uniform limit the HNC energy simplifies as: 

h 2 k 2 (S(k) - 1) " 



EuM 1 9 + 1 ( 



N 2 v ' 2 J (2vr) 3 p 



(S(k) - l)V(k) + 



(33) 



4m S(k) 

an expression that can be readily obtained from Eq. by assuming g(r) ~ 1, so that 
In g(r) w g(r) — 1. The last term in the integral is the kinetic energy contribution, and by 
comparing it with (|25|) one readily realizes that the momentum distribution is related to the 
static structure function in the UL through the relation: 

The condensate fraction, no, can then be recovered by imposing the normalization condition 
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Minimization of Ejjl with respect to S(k) gives the Euler-Lagrange equation in the UL, 

which can be solved for the the static structure function to obtain: 

S UL (k) = m = . (36) 

^jt 2 {k) + 2t{k)V{k) 

This expression is formally identical to the HNC one of Eq. (|19)L where V p h(k) has been 
approximated by the Fourier transform of the bare potential V(r). 

We notice that Sjji{k) coincides with the static sturcture function given by the RPA 
approximation to the dynamic susceptibility, x{k, ^>)*^> 

XRPA(k,u) = , (37) 

1 - V(k)Xo[k,w) 

where xo{k,u) is the dynamic susceptibility of the free Bose gas at zero temperature, 

Xo(k,u) = } , } . . (38) 

to — t(k) + ie u + t{k) + xe 

1 /2 

The poles of XRPA(k,uo), ^RPA(k) = (t(k) 2 + 2t(k)V(k)j , define the excitation energies of 
the system, while the imaginary part of XRPA{k, ui) gives the dynamic structure function, 

1 t(k) 

S RPA (k, u) = ^SXRPA(k, u) = 7T\ 5 ( U - £RPA(k)) . (39) 

7T e RPA {k) 

Integration of Snp A (k,uj) over u> leads to the static structure function, SppA^k), 

S RPA (k) = m = = , (40) 

^t 2 {k) + 2t{k)V{k) 

which coincides with Sub{k). Moreover, the UL energy ()33|) can be obtained by adding 
the RPA correction, AE R pa, to the uncorrelated energy (|3~Tj). AErpa can be evaluated by 
performing a coupling constant integration, 

S\(k) being the RPA static structure function corresponding to the interaction XV (k). After 
a straigthforward integration the uniform limit expression is recovered, 

Eul = \v{$) + AE RPA . (42) 
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A similar result has been found in Ref.— , where the soft-sphere gas has been studied in 
the framework of the HNC theory with an independent pair correlated wave function. The 
authors have shown that neglecting the composite diagrams and summing the chain diagrams 
only, the HNC/IPC approach leads to the RPA (and UL) energy functional (|33j) . All three 
theories produce the same description of the soft-sphere gas in the low density regime, since 
both uniform limit and correlated theories, in conjuction with a proper minimization via 
the Euler equation, take into account the long range correlations relevant in this region and 
correctly considered by the RPA. 



A. Some aspects of the low density expansion 

The formalism outlined in Section II allows to describe the ground state of any boson 
liquid or gas characterized by a hamiltonian of the form given in Q. At very low densities, 
however, all interacting Bose gases having the same scattering length are expected to behave 
in a similar way, as originally pointed out by Bogoliubov'. The energy of a Bose gas follows 
the universal form ((!]) as long as the system is dilute. However, the overall energy comes 
from a balance between kinetic and potential contributions. For hard spheres the energy 
is entirely kinetic, whereas for soft spheres the relative contribution of the potential energy 
may be as large as possible, according to the softness of the interaction. 

Despite the universality exhibited by the total energy, it is not clear that other ground 
state properties may show an analogous behavior. An answer to this question can be ob- 
tained by analyzing the hard and soft sphere gases and comparing with the predictions 
obtained within the Bogoliubov model. For instance, from the Bogoliubov excitation spec- 
trum, 



e B (k) = \Jt 2 (k) + 2t(k)4:7iah 2 p/m , (43) 
an approximate static structure factor can be obtained assuming a Feynman-like spectrum, 

S Bik) = Jg . (44) 

However, g(r) extracted from Ssik) is unphysically divergent at r = 0. 

In the Bogoliubov model the momentum distribution, rig (A;), at k > can be written in 
the forirA 

A 2 

M*)=i3i2 > Ak = 1 + ^r P m ~ €B{k)] ■ (45) 
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This expression coincides with n(k) as obtained in the UL (Eq. (|34j) ) when V(k) is replaced 
by Anaph 2 /m. 

Normalization of the full momentum distribution gives the fraction of particles in the 

k = state, n s, 

"OB = 1 " ■ ( 46 ) 

The kinetic energy computed through ne(A;) is divergent since nsik) ~ k~ A at large k. 



IV. DILUTE HARD SPHERES 



In this section variational results for the energy, two-body distribution function, static 
structure function and one-body density matrix (and momentum distribution) of a dilute gas 
of hard spheres are shown and discussed. The driving quantity of any variational approach 
is the total energy, given by the sum of the kinetic and potential terms. Since inside the core 
of the potential g(r) = 0, the energy is purely kinetic and, in units of h 2 /2ma 2 , becomes 

E(x) (f) 



~x J drg{r)V 2 \nf{r) . (47) 



N N 

We will adopt the long ranged correlation function obtained by solving the Euler- 
Lagrange equations and will analyze the related asymptotic behaviors. However, the HNC 
results for a simpler correlation function having a short-range structure, fsii{ r ), are also 
discussed. fsii(r) is choosen in such a way to minimize the lowest-order in the cluster ex- 
pansion of the energy of the homogeneous gas of HS with a healing condition at a distance 
d : taken CIS db variational parameter—. In the HS case fsRir^) is. 

f r < 1 

fsn(r) = {„.,„, ni (48) 

d sm[K(r-l)] i 
^ r sm[K(d-l)} ' ^ L ' 

where distances are in units of a and K fulfills the equation: cot[K(d — 1)] = (Kd)^ 1 . The 
latter condition ensures the healing properties: fsn(r > 1) = 1 and f' SR (r = d) = 0. 

The scaled energies per particle, E(x)/N, of the HS gas as a function of the gas parameter, 
x, are given in Table HJ Eel is obtained by solving the Euler-Lagrange equation, while in 
computing Esr the short-range correlation (J48|) has been used. Both energies are computed 
within the HNC/0 approximation, justified a priori by the smallness of x, and a posteriori 
by the eventual agreement with the exact DMC results 1 , reported also in the Table together 
with the results of the LD expansion ((TJ). 
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Eel and Esr are upper bounds to the exact DMC energy. Furthermore, since the solu- 
tion of the Euler-Lagrange equations yields the minimum energy for a Jastrow type wave 
function, the inequality E DMC < E EL < E S r holds at each density. Violations of this 
hiearchy at low x-values are probably due to numerical inaccuracies, rather than to the 
HNC/0 approximation. 

The low-density expansion does not satisfy the upper bound property. However, the 
lowest order of the same expansion, Eldq/N = 4irx, is a rigorous lower bound to the exact 
energy^. The overall differences among the EL, SR and DMC energies are very small (at 
most 6% in the worst case at x = 0.1) and both the EL and the SR results can be taken as 
good estimates. Using the HNC/0 scheme is therefore well justified in the range of densities 
explored, especially at the lowest ones. For the sake of comparison, a Variational Monte 
Carlo (VMC) calculation at x = 0.1, with fsR{r) having a healing distance d = 6, has 
also been carried out in order to estimate the relevance of the elementary diagrams. The 
result E V Mc/N = 3.74 ± 0.02 is to be compared with E SR /N = 3.97 in HNC/0. Hence, 
at the highest density considered the elementary diagrams contribute by only ~ 6% to the 
energy. Further differences with the DMC results have to be attributed to deficiencies in 
the two-body correlation factor and to the lack of three- and higher body correlations in 
the trial wave function. The energies in these approximations are plotted in Fig. ^ where 
the subtle differences between the points are hardly appreciable, except at the highest x's. 

The influence of the optimization on the energy is rather small. The energy is dominated 
by the short range structure of the potential, which requires the two body distribution 
function to be zero inside the hard core. The Euler-Lagrange procedure is instead important 
in establishing the long range structure of the distribution function g(r), or, alternativeley, 
of the low k behavior of S(k). 

The EL radial distribution function is shown in Fig. El for several values of the gas pa- 
rameter. At low x's g(r) is a monotonically increasing function of the distance, approaching 
faster and faster the asymptotic limit, g(r) — > 1, with the density. This behavior is readily 
understood recalling that g(r) measures the probability of finding two particles at a distance 
r, and the average interatomic spacing decreases with the density. At large densities the 
radial distribution function develops a local maximum close to the core radius. 

The radial distribution functions, gEii^), solutions of the EL equation at x — 0.001 and 
x = 0.005, are compared in Fig.|3]with the corresponding short-range ones, gsR( r ), computed 
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from fsn{r). At these densities, the energy of the system is accurately described by the 
expansion of Eq. (fl|). and it is natural to compare other ground state quantities to those 
corresponding to the Bogoliubov approximation. Fig. |H] also shows the radial distribution 
function, ^(r), obtained as the Fourier transform of Ssik), 

~ Vk 4 + 167rf = ' ^ 
QELij") and gsnij) are close at short distances: they both vanish inside the core and 
approach similarly the unity. However, differences with Qb{^) are significant. Qb{t) becomes 
unphysically negative at short distances to finally diverge at r — > 0. In fact, Ss{k) realisti- 
cally reproduces only the low k behavior of the static structure function in dilute systems. 
This corresponds to the large-r region in the associated radial distribution function. For the 
same reason, g_e(r) never develops a maximum. The three radial distribution functions have 
different asymptotic behaviors, as shown in Fig. 0J where r 4 (g(r) — 1) is given in the EL, SR 
and Bogoliubov approaches at x = 0.005. Both gEL{ r ) — 1 and gB( r ) — 1 behave as r -4 at 
large r, a property not shown by gsR(.r) — 1 since fsn( r ) does not have the appropriate long 
range behavior, f(r — > oo) — 1 ~ r -2 . Actually, the long range behavior of g(r) is related to 
the sound velocity, c (in units of h/2ma), by 9 

11 

g(r -> oo) ~ 1 , (50) 

n 2 cx r 4 

and consequently, the static sturcture function goes like 

S(k^0)^-. (51) 
c 

In the Bogoliubov approximation, 

1 1 

g B {r -> oo) ~ 1 - A7{5/2xV2 - A , (52) 



and cb = vT&nx, consistent with the value provided by the compressibility, k t , c = 
(pttikt) -1 ~' 2 / (h/2ma), by keeping only the first term in the low-density expansion of the 
energy. At x = 0.005 the constant value of r A (gB{r) — 1) = —40.4 gives cb = 0.50, smaller 
than the estimated cel = 0.61. 

The short-range radial distribution function is compared in Fig. 03 with the VMC one 
at x = 0.01 and 0.1, for the same two-body correlation factor. We notice two aspects: 
(i) the inclusion of the elementary diagrams enhances the peak of the radial distribution 
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function, and (ii) there is a remarkable difference between fsn{r) and g(r). The many-body 
contributions included by the HNC scheme cannot be neglected and, therefore, g(r) cannot 
be approximated by its lowest order cluster expansion value, gLo{ r ) — f 2 { r )- 

The optimal structure function S(k) is shown for several values of x in Fig. El The slope 
of S(k) becomes smaller when x increases because the speed of sound increases with the 
density. Furthermore, the linear low k behavior, which guarantees the correct low energy 
excitation spectrum, is evident. However, the linearity of the static structure holds only at 
very low k, as can be seen from Fig. [7J where the ratio S(k)/k is shown in both the EL and 
the Bogoliubov cases at x = 0.01. In the later case the linear regime is valid only when 
k < VT6nx = c B (see Eq.flHl)). 

The particle-hole interaction V p h (Eq. (J2H|) ). is shown in Fig.|H]for x = 0.001, x = 0.05 and 
x = 0.1. The left panel gives V p h in r-space. Two regions are separated by a discontinuity at 
r — I, produced by the term | V 'J 1 g{r) | 2 in Eq. (J2U|) since the first derivative of the HS radial 
distribution function is discontinuous at the core. As already noticed, — V P h(r < 1) coincides 
with the induced interaction, uj(r) ()21)1 . At the lowest densities | V \J ' g{r) | 2 dominates 
and the strength of V p h{r) is almost entirely exhausted by it. The right panel displays the 
Fourier transform of V p h at the same densities. V p h(k) is an oscillating function of k, with its 
highest amplitude at k — 0, and V p h(k = 0) = c 2 /2. Notice that the figure shows V P h(k)/x, 
therefore V p h } B(k)/x would be 8it constant and independent of k and x, which turns out to 
be a very bad approximation to V p h{k)/x. However calculating the speed of sound from the 
two terms of Eq. (0), one gets c 2 /2x = 8tt + Yl&^^x 1 ! 2 ^ in much better agreement with EL 
V ph (k = 0)/x. 

The Bogoliubov estimate of the particle-hole interaction fEqs.(|T^|) and (jinj)) is V p ^B{k) = 
8ttx and it is a rather poor approximation to V p h{k). In this case, the excitation spectrum, 
e(k) = t(k)/S(k) = \Jk A + 2k 2 V p h{k), does never develop a rotonic structure, since V p h,B{k) 
is a constant. In the EL approach, instead, the first oscillations of V p h{k) (and of the static 
structure function) are large enough to produce a maxon-roton behavior at the largest x 
value (x = 0.1). The two spectra are explicitely shown in Fig. El The Figure also gives e(k) 
at x = 0.001, where the differences between the EL and Bogoliubov results can hardly be 
resolved. 

The last quantity analyzed is the momentum distribution, n(k). Table HP reports the con- 
densate fraction, n , and the kinetic energy for both the EL and the SR correlation factors. 
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The Table gives the kinetic energies estimated by Eq. (J25j) . T n , and by the second term of 
Eq. (JT5j) . T 9 . At low x values the two estimates almost coincide, whereas, discrepancies at 
larger densities are to be ascribed to the lack of the elementary diagrams contribution in 
the cluster expansions of the radial distribution function and the momentum distribution. 
It has to be noticed that these diagrams play different roles in the two expansions^. 

The condensate fraction decreases with x, running from ~ 95% at x = 0.001 to ~ 33% 
for x = 0.08. The corresponding predictions for ti b (Eqflfij) are n s(x = 0.001) = 0.95 and 
Uqb{% = 0.08) = 0.57. This points again to the failure of the Bogoliubov model at large 
densities. 

The momentum distributions in the EL, SR and Bogoliubov cases at x — 0.05 are shown 
in Fig. The optimal n(k) has the long-wavelength limh-& 

\im kn(k) = , (53) 
k^o 4 

while ns{k) satisfies an analogous relation, with uq = 1 and c = cb- In contrast, nsR.(k) 
does not behave as 1/k at k — > 0, and kn(k) vanishes at the origin. This fact is due to the 
lack of the proper long range structure in f$R- 

Fig. ^2 shows the optimal momentum distribution at x = 0.08 and 0.05. At the lowest 
x values kn(k) is a monotonically decreasing function of k. Actually, this is also the case 
in the Bogoliubov approximation at any density. When x increases, kn(k) develops a peak 
in the EL case at low k. This maximum can be considered as a genuine effect of the short 
range correlations in the EL approach. Notice that the value at the origin of kn(k) results 
from a competition between uq and c. These two quantities largely vary with the density (n 
decreases and c increases with increasing density), but the overall variation of linn^o kn(k) 
is less pronounced. 

V. DILUTE SOFT SPHERES 

The soft spheres potential (Eq. ©) is characterized by a core radius R and a potential 
height Vo that determine the scattering length a according to Eq. (jlOj) . The SS scattering 
length is always smaller than the core radius, approaching it when Vq increases. In the 
very dilute regime, HS and SS systems, having the same scattering length, are expected to 
have the same energy, while their separate kinetic and potential contributions may differ. 
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Moreover, at low x-values the total energy is well reproduced by the two terms of the low 
density expansion of Eq. (0). In this section, we will study the deviations from this behavior 
and the influence of the shape of the potential on the energy and its components. To this 
aim, two SS potentials with the same scattering length and different radii are considered: 
the SS10 and SS5 potentials, having R = 10 and R = 5 in units of a, respectively The 
corresponding heights are Vq SW = 0.00681670 and Vq S5 = 0.06308561 in units of h 2 /2ma 2 . 

Table IIHI reports the scaled total, kinetic and potential energies per particle for SS10 
and SS5 in the EL, SR and UL cases, compared with the available exact DMC results^ at 
x = 0.0001 and x = 0.01. The low-density expansion yields E LD /N(x = 0.0001) = 0.001317 
and E LD /N(x = 0.01) = 0.1862. The lower bound energies are E lh {x = 0.0001) = 0.001257 
and Ei b {x = 0.01) = 0.1257. The upper bound energies (J3TJ) depend on the shape of the 
potential, and, for SS10 and SS10, give: E^ b sw /N(x = 0.0001) = 0.0014277, E^ 5 /N(x = 
0.0001) = 0.0016516, E^ 10 /N(x = 0.01) = 0.14277, and E^/N(x = 0.01) = 0.16516. 
Both, the upper and the lower bound properties are here satisfied. The shape dependence 
is weaker at low x, in fact the difference between the S'S'10 and SS5 energies is smaller 
at x = 0.0001 than at x = 0.01. On the other hand, the kinetic and potential energies 
exhibit a potential dependence at any x-value. The IPC results of Ref.— are also given 
in the Table. The EL energies are very close to the DMC ones at x = 0.0001 and the 
variational hiearchy is always fulfilled (DMC<EL<SR). As expected, the UL works better 
for the weaker SS10 potential, as well as the IPC wave function. The Jastrow wave function 
is variationally preferred to the IPC one at any x, as shown by the comparison between 
the EL and IPC energies. The good agreement between T g and T n seems to show that the 
HNC/0 approximation in the EL approach does not invalidate this conclusion. 

The behavior of the scaled energy per particle along x is shown in Fig. in units of 47rx, 
for the SS10, SS5 and HS potentials, in the EL, IPC and DMC cases, and compared with the 
upper bounds. The lower bound equals to unity in this units. Starting at x > 10~ 3 , clearly 
appears a shape dependence of the energy with the harder SS5 potential energies closer to 
the HS ones. Moreover, it is also to be stressed a dependence on the quality of the wave 
function in the same region, since differences between the Jastrow and IPC cases become 
evident for SS5. The quality of the upper bound improves when R increases, because V 
decreases at fixed a and the perturbative expansion is expected to converge faster. 

The energy and, in general, the structure of the ground state depends, to a large extent, 
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on the two-body correlation factor employed. In the EL and UL cases, f(r) is a derived 
quantity as the optimization is accomplished by varying g(r) or S(k). In the SR case, /5_r(t) 
is an input function containing some variational parameters. In the present work for the SS 
gas we employ 



which is flexible enough to obtain a reasonable value of the energy. Notice that the SS 
potential does not require the correlation function to be zero at the origin, so b ^ 1. 

Figure Uni shows g{r) for SS5 at x — 0.001 in the EL, UL and SR cases and /J^t"). 
Many-body effects make gEL^) quite different from both at low and intermediate 

r-values. The EL radial distribution function is softer at the origin, and, in general, less 
repulsive than the short-range one. The radial distribution function in the uniform limit 
approximation is close to f^ L at small distances, approaching gEL{f) at r ~ 5. The radial 
distribution function obtained by the Bogoliubov approach (not shown in the figure) would 
exhibit, as for the HS case, an unphysically divergent behavior at short distances. 

The dependence of the radial distribution function on the shape of the potential and on 
x is illustrated in Fig. El containing the EL distribution functions for SS5 and SS10 at 
x = 10~ 4 and x = 10 -3 . The central hole of the radial distribution function is obviously 
deeper for the more repulsive SS5 potential. It becomes less pronounced at higher x, since 
the probability of finding two bosons at short relative distances increases with the density. 
The universal behavior in x is recovered at large distances, where the Bogoliubov approach 
becomes reliable. 

The static structure functions, in the EL, UL and Bogoliubov approaches, are plotted in 
Fig. EHfor the SS5 interaction. In all cases, S(k) grows linearly at the origin, and the slope, 
governed by the sound velocity (Eq. (fBTjl ). is similar in all three cases. The Bogoliubov 
sound velocity, cb = \flQixx, is smaller than the UL estimate, 



and consequently, the slope at low momenta in Ss(k) is slightly larger than in Sui(k). 
Suhik) approaches SEh{k) and the asymptotic value faster than Sb(/c). Actually, at a given 
density Ss{k) and cb are identical for all the SS interactions with the same scattering length. 
The short-range g(r) does not have the proper asymptotic behavior and produces a static 
structure function that does not vanish at the origin and does not increase linearly at low 



fsR(r) = 1 - be 



—cr 



.2 



(54) 




(55) 
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k. S EL {k) at x = 0.0001 and x = 0.001 are shown in Fig. M for SS5 and SS10. The SS10 
sound velocity is smaller and the slope of S(k) at the origin is larger. The differences are 
more evident at the largest x. 

The particle-hole interaction corresponding to SS5 at x — 0.001 is shown in Fig. \H\ 
for the EL, UL and Bogoliubov approaches. The left and right panels give V p h(r) and its 
Fourier transform, V p h{k), respectively. VR(r) = 8irx5(r) is not shown in the figure. The 
p/i-interaction in the UL limit coincides with the interaction itself, Vp^ L (r) = V(r). V p E h L (r) 
is discontinuous at r = R. The differences between V p u h L {k) and V p EL (k) are larger at low 
momenta, consistently with the differences observed in the static structure function. 

Finally we discuss the results for the momentum distribution, n(k). The quantity kn(k) 
at x = 0.001 is plotted in Fig. [181 for the SS5 interaction in the EL, Bogoliubov, UL and 
SR cases. All momentum distributions but the short-range one satisfy the long-wavelength 
behavior (J33J). Since n SR (k = 0) ^ 0, knsR^k) = at the origin. As in the HS case, 
knsn{k) presents a maximum at intermediate momenta, as a byproduct of the absence of a 
long-range structure in the two-body correlation factor fsii(r). 

In the x-range considered, knEi{k) is a monotonously decreasing function of k. At higher 
x it develops a maximum at low momenta, consistently with the HS results, also found in 
the high density homogenous atomic 4 He. In the UL, 



and the value at the origin depends on x through the sound velocity, whereas its slope, 
—1/2, is independent on the density. 

kriB(k) is also a decreasing function of k, its k — value being slightly smaller than the 
UL one, corresponding to a lower sound velocity. The UL momentum distribution for the 
SS potential decays as k~ 8 , and the condensate fraction can be computed by imposing the 
normalization condition. The condensate fractions and the kinetic energies obtained from 
the momentum distribution are reported in Table ITTT1 

The dependence of n{k) on the shape of the potential is studied in Fig-HHJ where kriEL{k) 
is plotted at x = 0.0001 and x = 0.001 for SS5 and SS10. Apparently, the SS5 and the SS10 
momentum distributions are identical at x = 0.0001, pointing to a dependence only on the 
scattering length. This conclusion must be taken with caution, since in Table ITTT1 more than 
a factor of two between the SS5 and SS10 kinetic energies, T n , is found even at this low 
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density. 

The condensate fraction also depends on the shape of the potential. Figure 1201 shows 
n EL and tiq^l at x — 0.001 as a function of the radius of the SS potential, R, at a fixed 
value of the scattering length. As expected, the condensate fraction grows with R, since the 
interaction softens. The UL approach becomes accurate at large i?-values, in the very weak 
interaction limit. In all cases, no,ui is smaller than uo,el- 

VI. SUMMARY AND CONCLUSIONS 

We have carefully investigated the energy and structure of a homogeneous gas of bosons 
interacting via hard and soft sphere potentials. We have adopted and compared several 
many-body approaches, as the variational correlated theory, the Bogoliubov model and the 
uniform limit approximation, valid in the weak interaction regime. When possible, the 
results have been compared with the exact diffusion Monte Carlo ones. A Jastrow type 
correlation appears to produce a good quality wave function if the hypernetted chain energy 
functional is freely minimized and the resulting Euler equation is solved. This is true for 
both hard and soft sphere interactions. The study of soft sphere potentials has confirmed the 
appearence of a shape dependence in the energy at x ~ 0.001, as already found by the IPC 
calculations of Ref.— and by the DMC study of Ref.— . We have numerically compared the 
EL results with those obtained in the uniform limit of weak interaction, in the Bogoliubov 
approximation and in the IPC theory. As expected, the differences are more relevant for the 
strongest interactions. The uniform-limit and independent-pair-correlation energies become 
more and more reliable as the interaction weakens. The universality breaks at much lower 
a;-values for quantities different from the energy. For instance, the short range structure of 
the radial distribution functions (and, in consequence, the large momentum behavior of the 
static structure function) largely depends on the shape of the potential already at x — 0.0001. 
We find a potential shape dependence in the SS condensate fraction at x — 0.001. These 
results, as well as those for the energy, may help in evaluating the parameters entering the 
nonuniversal corrections in effective field theory^, whose extraction from the available DMC 
calculations is heavily biased by the statistical errors. 

The presence of a maximum in the correlated distribution function when x increases has 
been discussed in details. The Bogoliubov model does not show this short range structure, 
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independently on x. The excitation spectrum develops a rotonic maximum at large x, related 
to the shape of the particle-hole interaction. Again, the Bogoliubov approach, providing a 
constant V p h{k), does not allow for a roton-like excitation. The effect of the correlations 
is also apparent in the momentum distributions. In the hard sphere gas the CBF kn(k) 
acquires a maximum at low momentum at x ~ 0.06. Both the Bogoliubov and the uniform 
limit (for soft spheres) approaches fail to reproduce this behavior. 

The correlated basis functions theory may be extended to treat fermionic hard and soft 
spheres or mixtures of Fermi-Bose gases. Work along these lines is in progress. More chal- 
lenging is the application of the CBF approach to trapped atomic gases, without resorting to 
local density type approximations^. Employing the variational method in finite dilute sys- 
tems presents mainly technical difficulties. However, it will probably represent the natural 
development of this technique. 
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FIG. 1: Scaled energy per particle of the HS gas as a function of x. Solid circles: Euler-Lagrange 
results; stars: HNC/0 results for the parametrized two-body correlation factor of Ea.(|48|): open 
squares: Diffusion Monte Carlo resuls; dashed line: first two terms of the low density expansion, 
Eq.©. 
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FIG. 2: EL radial distribution function g(r) for HS at x = 10 4 (short-dashed line), x = 10 3 
(long-dashed line), x = 10 -2 (dot-dashed line) and x = 10 _1 (solid line). Distances in units of a. 
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FIG. 4: Long range structure of the HS g(r) at x = 0.005 in the EL (solid line), SR (dashed line) 
cases compared with the r — > oo contribution in the Bogoliubov approximation (dot-dashed line). 
Distances in units of a. 
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FIG. 5: gvMc( r ) (solid line) and gHNC/o( r ) (dashed line) computed from fg R (r) of Eq.(@HJ (dot- 
dashed line) for HS at x = 0.01 (left panel) and x = 0.1 (right panel). Distances in units of 
a. 
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FIG. 6: HS EL static structure factor S(k) at x = 10 4 (solid line), x = 10 3 (dotted line), 
x = 5 • 10~ 3 (dashed line) and x = 10~ 2 (dot-dashed line). Momenta in units of a -1 . 
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FIG. 12: Scaled energy per particle for the SS5 (triangles) and SS10 (circles) potentials in the EL 
(filled triangles and circles) and IPC (empty triangles and circles). The stars and crosses are the 
available DMC energies for SS10 and SS5, respectively. The solid line represents the EL energies 
for the HS potential. The horizontal lines give the upper bounds for the SS5 (dash-double dotted) 
and SS10 (dashed) potentials. In these units, the common lower bound equals 1. 
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FIG. 13: Two-body distribution functions for the SS5 potential at x = 0.001. The solid, dot- 
dashed and dashed lines correspond to the EL, UL and SR cases, respectively. The two-body 
correlation function, / 2 (r), as extracted from the EL g(r) is plotted as a dotted line. Distances in 
units of a. 



36 



WO 




FIG. 14: EL radial distribution functions for SS10 and SS5 at x = 0.0001 (solid line for SS10 
and dotted line for SS5) and x = 0.001 (dashed line for SS10 and dashed-dot-dot line for SS5). 
Distances in units of a. 
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FIG. 16: EL static structure functions for SS10 and SS5 at x = 0.0001 ( solid line for SS5 and 
dot-dashed line for SS10), and x = 0.001 (dashed line for SS5 and dotted line for SS10). Momenta 
in units of a -1 . 
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FIG. 17: Particle-hole interaction in r-space (left panel) and in k-space (right panel), for SS5 at 
x = 0.001. The solid and dot-dashed lines correspond to the El and UL cases, respectively. The 
dashed line in the right panel gives the Bogoliubov V p h- Distances in units of a and momenta in 
units of a -1 . 
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FIG. 18: kn(k) for SS5 at x = 0.001. The solid, dashed, dash-dotted and dotted lines correpond 
to the EL, Bogoliubov, UL and SR results respectively. Momenta in units of a -1 . 
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FIG. 19: EL kn(k) at x = 0.0001 (solid line for SS10, dashed line for SS5) and x = 0.001 (dotted 
line for SS10, dot-dashed line for SS5). Momenta in units of a -1 . 
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FIG. 20: Condensate fraction for the soft sphere gas, at x = 0.001, as a function of the radius, 
R, of the SS potential at fixed scattering length. EL, solid line; UL, dashed line. The Bogoliubov 
condensate fraction is uqb = 0.952. Distances in units of a. 
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5-10" 2 1.338 1.383 1.402 7 1.305 

0.1 3.627 3.848 3.971 6 3.170 



TABLE I: Scaled energy per particle for the hard spheres model, as a function of x. See text, d 
is the healing distance corresponding to the correlation factor of Ea.(|48[). in units of the scattering 
length. In Eld the first two terms of the expansion ^ have been used. 
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TABLE II: Scaled kinetic energy per particle of the HS gas calculated using the momentum dis- 
tribution and the radial distribution function for the EL and SR wave functions at several values 
of x. no, n-QDMc an d uqb are the EL, DMC and Bogoliubov condensate fractions. 
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-4 


0.992 


IPC 


io- 


-4 


10 1.311 • 


10" 


-3 






















DMC 10" 


-4 


10 1.303 • 


10" 


-3 




















0.989 


EL 


10" 


-1 


5 1.314- 


10" 


-3 


1.044 ■ 


10" 


-3 


2.708 • 


10" 


-4 


2.630 • 


10" 


-4 


0.985 


SR 


10" 


-1 


5 1.361 • 


10" 


-3 


1.138 ■ 


10" 


-3 


2.231 • 


10" 


-4 


2.231 • 


io- 


-4 


0.996 
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UL 


10" 


-4 


l.Zol • 


1 n- 

1U 


-3 


U.oDo • 


i ri- 
ll) 


-3 






O. f oU 


1 n~ 
• ill 


-4 


u.yoz 


IPC 


10" 


-4 


l.ool • 


1 n- 


-3 


















n 087 

u.yo t 


DMC 10- 


-4 


o l.ouy 


1 n- 
1U 


-3 


















u.yoy 


EL 


io- 


-2 


10 1.404 • 


1 n- 
10 


-1 


1.394 • 


10 


-1 


n nnn 

O.yyO • 


10 


-3 n noi 

o.ysi 


1 n- 
10 


-3 


o.yso 


fin 

bR 


10 


-2 


1U 1.4UD • 


1 n- 

1U 


-1 


i.oyo • 


1U 


-1 


u.you • 


1 ri- 


-3 n QRn 
u. you 


i ri- 
ll) 


-3 


n 0.77 

u.y < / 


UL 


10" 


-2 


in i /in/i 

1U 1.4U4 • 


1 n- 

1U 


-1 


i.oyo • 


i ri- 
ll) 


-1 






n qrq 
u.yoo 


1 ri- 
ll) 


-3 


n osn 
u. you 


IPC 


10" 


-2 


1U 1.4Uo • 


1 n- 

1U 


-1 




















EL 


10" 


-2 


5 1.532 • 


10" 


-1 


1.468 • 


io- 


-1 


6.445 • 


10" 


3 6.480 


10" 


-3 


0.951 


SR 


10" 


-2 


5 1.535' 


10" 


-1 


1.481 • 


io- 


-1 


5.350 • 


io- 


3 5.350 


10" 


-3 


0.960 


UL 


10" 


-2 


5 1.528' 


10" 


-1 


1.464 • 


io- 


-1 






6.375 


10" 


-3 


0.950 


IPC 


10" 


-2 


5 1.556' 


10" 


-1 


















0.953 



TABLE III: Scaled energies and condensate fractions for the soft sphere gas in the El, SR, UL, IPC 
and DMC approaches for two potentials having the same scattering length (see text) at x = 10~ 4 
and x = 10~ 2 . 
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